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Abstract 

In both nature and engineering, loosely packed granular materials are 
often compacted inside confined geometries. Here, we explore such be¬ 
haviour in a quasi-two dimensional geometry, where parallel rigid walls 
provide the confinement. We use the discrete element method to inves¬ 
tigate the stress distribution developed within the granular packing as a 
result of compaction due to the displacement of a rigid piston. We observe 
that the stress within the packing increases exponentially with the length 
of accumulated grains, and show an extension to current analytic models 
which fits the measured stress. The micromechanical behaviour is studied 
for a range of system parameters, and the limitations of existing analytic 
models are described. In particular, we show the smallest sized systems 
which can be treated using existing models. Additionally, the effects of 
increasing piston rate, and variations of the initial packing fraction, are 
described. 


1 Introduction 

When granular materials are placed in confined geometries, we often observe a 
significant portion of the stress being redirected towards the confining bound¬ 
aries. This phenomenon has been systematically studied for many systems 
miniiis], most notably in silos, beginning with [5]. Force redirection has been 
attributed to the granular nature of the material, and has in many cases been 
shown to be well represented by a constant coefficient, known as the Janssen 
coefficient K, defined in one spatial dimension as 

K = 

where (Jr is the redirected stress due to some applied normal stress cr„. Here 
we investigate the development of stresses within a granular packing, confined 
between two horizontal plates, subjected to a rigid piston impacting it from 
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one side. As the piston moves, granular material is compacted near the pis¬ 
ton, and with increasing displacement of the piston, the size of the packing 
increases. Such an accumulation process is known to occur in the petroleum 
industry, where sand is liberated from the host rock during extraction, altering 
the underground morphology of cracks [iHiiin]. This may also be relevant for 
understanding proppant flowback in propped fractures m- Additionally, this 
geometry is representative of a number of recent experimental studies in Hele- 
Shaw cells m El [HITS] where the validity of Janssen stress redirection has not 
been ascertained. 

There are a number of interesting patterns which form when a granular 
material is displaced by a flexible interface in such a geometry El [15]. The 
nature of the patterns have been shown to depend on many factors, primarily 
the initial packing fraction and the rate of displacement [14]. For this reason, we 
here investigate the microstructural and mechanical evolution of such a system 
under these conditions. To reduce the complexity of the system, we consider 
only a rigid piston. 

We are interested in systems which are highly confined. In common experi¬ 
ments with granular material inside Hele-Shaw cells, there are in general fewer 
than 20 grain diameters between the two Hele-Shaw plates, typically down to 
around 5 grain diameters El- As the confinement increases, i.e. as the gap 
spacing decreases, we expect a transition from three dimensional behaviour to¬ 
wards a behaviour governed by the boundaries, as demonstrated for vertical 
silos in [T] . It is then of interest to study the changes that result from increas¬ 
ing confinement. We expect that altering the confinement will affect the force 
redistribution. A transition may occur for extremely confined systems where 
such an assumption concerning force redirection may not be valid. 

Existing analytic models for the micromechanics of such a system generally 
reduce the problem to one spatial dimension (a;), assuming that variations in 
both remaining directions {y and z) are small, although recently curved inter¬ 
faces have also been described [4]. For simplicity, we consider a flat interface, 
and validate the analytic description with a discrete element model. 

Firstly, in Section describe the numerical model that has been used to 
simulate this system. In Section [^ we establish continuum properties which 
correspond to the analytic formulation, and show comparisons between the two. 
In particular, the limitations of current analytic models are identified. Finally, 
a parameter set is proposed that best fits the analytic theories for a wide range 
of system variables. 

2 Numerical model 

This paper is an investigation into the micromechanics of a system which is 
highly constrained by external boundaries. For this reason, it is ideal to use a 
particle based method to model the behaviour, as the total number of grains 
in the system is small. Towards this end, we use a conventional soft sphere 
discrete particle approach, implemented in the open source code MercuryDPM 
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Figure 1. Particle positions during a single test, shown at six piston dis¬ 
placements s. Top to bottom: s = 0,10,20,30,40 and 50. Labels refer to the 
coordinate system x, y and z, the gap height £>, the plug length L and the initial 
packing fraction 4>o- Particles are coloured by size, darker colours representing 
smaller particles. 


(www.mercurydpm.org) [niEi]. The geometry considered here, shown in Figure 
consists of a rigid piston, oriented in a space r = {a;, y, z}, with normal along 
the X axis, which pushes particles between two rigid walls, separated by a spacing 
D and having normals in the ±z directions, with two periodic boundaries in the 
remaining perpendicular direction, y. The coordinate system moves with the 
piston, such that it is located at a: = 0 at all times. As the piston moves 
horizontally at a velocity u towards the grains, its displacement at any time t 
is then s = ut. 

We work in a system of non-dimensionalised units with the following proper¬ 
ties; length and mass have been non-dimensionalised by the length and mass 
of the largest particle in the system, respectively, where the prime indicates 
that the quantity has dimension. The particle diameters, d, we use are therefore 
d < 1, with material density defined by 47r(l/2)^pp/3 = 1, or pp = Q/tt. Time is 
non-dimensionalised by the time taken for the largest particle to fall from rest its 
own radius under the action of gravity, so that a unit time is t = ^dm/g, which 
requires that g = 1. Other values are non-dimensionalised by a combination of 
these three scales, for example stress is non-dimensionalised by m^y'/d'^. 

Particles are filled into the available space by assigning them to positions 
on a regular hexagonally close packed lattice, dimensioned such that particles 
of diameter 1 would be in contact. In all cases we use particles of 0.5 < d < 1 
to avoid crystallisation. Variable particle filling is facilitated by changing the 
number of layers of the grid, such that the initial packing fraction, (pQ defined in 
Eq is approximately constant throughout the cell. From t = —10 to t = 0, 
gravity in the —z direction is increased from y = 0 to y = 1 to settle the grains 
in a loose packing. From t = 0 the piston begins to move at velocity u. 
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Figure 2. Contact laws used in the discrete element model. Interactions are 
characterised by (a) normal, (b) tangential and (c) rolling laws, parameterised by 
stiffnesses k, kt and k^, viscous dissipation 7 , 7 ^ and 7 ^ and friction coefficients 
Ht and ^r- Full details are given in El- 


Table 1. Material properties of the spheres. 


Direction 

Stiffness, k 

Dissipation, 7 

Coefficient of friction, p, 


( kg/s= \ 

( ^ 



yrnLKa'dD ) 



Normal 

100000 

1000 

- 

Tangential 

80000 

0 

0.4 

Rolling 

80000 

0 

10-3 

As shown 

in Figure 

the particles’ material properties are described 


normal, tangential and rolling damped spring sliders El Hi! with the properties 
contained in Table These values have previously been calibrated to mimic 
micron sized silica beads [5]. The walls are implemented such that they are 
rough; when a particle contacts a wall, the piston, or both, it is prohibited from 
rotating, i.e. = 1. Otherwise, the interaction properties are the same as 
between two particles, except that the walls and piston are of infinite mass. We 
therefore have a well defined macroscopic sliding friction of ^ = 0.4 that does 
not depend on the rate of loading. 

In the following Section we will firstly detail the important macroscopic 
quantities measured from a single simulation. We will then investigate the effect 
of three controlling parameters on the evolution of the system: the Hele-Shaw 
spacing Z), the initial packing fraction (j)Q and the velocity of the piston, u. The 
piston rate, however, is not a priori a governing quantity, so we choose to control 
the piston rate via the inertial number, I, which is defined as / = jdml\/Pjpp, 
which is the ratio of inertial to imposed stresses, where P is a typical pressure 
and 7 is a typical shear strain rate [ 6 ]. Taking 7 = u/P, and P = PppD, gives 
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3 Results 

As depicted in Figure three distinct regions exist inside the system. These 
are termed the plug zone, where particles are densely packed near to the piston 
{x < L), the undisturbed zone, far from the plug, and the transition zone, where 
the plug is accumulating. To define these regions systematically, we must first 
measure the solid fraction, iz = 14/Vt, which is a local measure of the ratio of 
the volume of solids to the total volume. The solid fraction is coarse grained in 
one spatial dimension, x, as 


’^{x,t) = .^'^ViVVix - Xi{t)), (1) 

i—1 

where N is the number of grains in the system, W is the width of the system, 
Vi is the volume of the i-th grain, Xi its centre of mass and W is the coarse 
graining function, chosen in this case to be a ID normalised Gaussian function 
mm- Such a coarse graining method allows the coarse graining width to be 
defined, so that either the macro- or micro-structure is visible. We choose to 
use a coarse graining width equal to the maximum particle diameter, such that 
small scale variations are minimised, and smooth continuum fields are obtained 
|20j . All other continuum quantities are defined using the same coarse graining 
technique, presented in mm- Using the definition Q , we denote the maximum 
solid fraction, Vm, as an average over the solid fraction close to the wall at some 
time when the transition zone is far from the piston as 

2 j-WD 

i'm = ^ v{x) dx, (2) 

5D J^D 

where in practice a numerical integration is done over the discrete coarse grained 
cells, and the limits of 5D and 101? are chosen arbitrarily. We define the nor¬ 
malised packing fraction, (j){x,t), as (f> = vjvrm and note that in the absence of 
volumetric expansion or dilation of the granular material, (j) is directly propor¬ 
tional to the height of the packing between the Hele-Shaw walls. The undis¬ 
turbed zone is that which is maintained at the initial packing fraction c/iq, which 
is defined at time t = 0 as 


2 rWD 

^0 = FU / (t){x,t = 0) dx. 

5D J 5 D 

To delineate the plug, transition and undisturbed zones, at each time t we use 
linear regression to find the best linear fit to the points in the range (p = ^0 + 0.1 
to (j) = 0.9, such that pD k, b — xtan0, for some value of h and slope angle 9. 
Five examples of this are shown in Figure |3]4. Using this best fit, we can define 
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Figure 3. Evolution of coarse grained properties of the system with increas¬ 
ing piston displacement, s, for D = 5, (j>o = 0.5 and I = 0.01. (A) Five 

examples of the normalised packing fraction of the particles, (j). The magenta, 
red, green, blue and black lines indicate displacements which correspond to 
L = 0, 12.5, 25, 37.5 and 50 respectively. The same colour scheme is used for 
each subsequent plot in this Figure unless stated otherwise. A filled circle indi¬ 
cates the measured value of L, and the cyan dashed line indicates the linear best 
fit measurement of the transition zone. (B) Evolution of the measured value of 
L with increasing displacement s drawn in black. The cyan line indicates the 
best linear fit to these points. (C) The coordination number, Z, as a function 
of normalised distance from the piston head, x/L. (D) Normal stress measured 
at the piston is shown in black. The cyan dashed line indicates the best fit 
estimate of Equation ([^. (E) Normal stress distributions within the packing. 
Solid lines indicate a^x-, dashed cfyy and dotted (j^^- (F) Apparent friction co¬ 
efficient measured within the packing, ^ = \<Jxz\/<Xzz- (G) Out of plane Janssen 
coefficient measured within the packing, Ky = ayyjuxx- (H) By plane vertical 
Janssen coefficient measured within the packing, = {cFzz — Ppi^gD/Z)/Uxx- 
(I) Absolute value of the x component of the eigenvector of the major principal 
stress. 
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the point at which the plug zone meets the transition zone as the intersection 
of the best fit regression line and ^!> = 1, such that L = {h — D)/ tan0, as shown 
in Figure [3]4. The value of L grows monotonically with piston displacement, as 
shown in Figure]^. Conservation of mass implies that on average, if there is no 
volumetric change in the packing, and no slip between the grains and the walls, 
this relationship can be expressed as 


L _ 

s l-4>o’ 


(3) 


which is shown as the dashed line in Figure [^. For all cases reported here, the 
value of 9 does not appear to vary with increasing plug length L. The point at 
which the undisturbed zone meets the transition zone can then be defined in a 
similar manner as above, by the intersection of the best fit regression line with 
(j) = (j)o. The coordination number, Z, (average number of contacts per particle), 
is fairly constant in the plug zone, (Figure]^), increasing with compaction, as 
rearrangement occurs. Additionally, at large values of L the stresses are high 
enough to cause significant overlap of the particles (up to 1%). In the transition 
zone, significant particle rearrangement lowers the coordination number. 

Coarse graining techniques in general cause measured fields to converge to¬ 
wards zero near boundaries [21]. While it is feasible to reconstruct these fields 
in general near individual boundaries, near the piston we have three distinct 
boundaries which all interact. To access stresses in this region, it is then prefer¬ 
able to directly measure the forces applied to the rigid boundaries of the system. 
Towards this end, we denote as the normal stress measured at the piston. 
This is shown as a function of piston displacement in Figure]^. The stresses 
measured from coarse graining within the packing are shown in Figure [^. In 
both cases, the stresses grow exponentially both with increasing L and decreas¬ 
ing X, as shown by the linearity in a semilog space in Figure [dp and E. 

An analytic expression to describe this stress evolution was first derived in 
|9], assuming that: (a) the stress redirection in the z direction is described by 
a constant where — Dppiymg/‘2 is the component 

of the vertical stress not due to gravity and (b) friction at the side walls is 
M = o’^y/o’zz- It can be shown using force balance in the x direction that under 
these conditions, if there is no net acceleration. 


d(T XX 
dx 


D 


^XX 




L). 


By further assuming that the stress at x 


L is a constant, i.e. aj- = <Jxx{x = 


Ppp^mg \ 2fiK^{L-x)/D _ ^Pp^m9 / a\ 

2K, ) 2K, ' ^ ’ 

Previously, the threshold stress ut has been modelled as either estimated 
from experimental data |4|, or as a constant by assuming sliding of a wedge 
of material p ng. Here, we choose to model the threshold stress ctt as a 
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Figure 4. Limit equilibrium of the transition zone. For the pile to be displaced 
by a force two forces must be overcome; Tb, the basal traction due to the 
weight Wq of the green region, and the x component of the internal sliding 
traction along the assumed failure surface denoted by the dashed line, due to 
the normal component of the weight above the failure surface, Wi, of the blue 
region. 


00 dependent quantity by considering limit equilibrium of a wedge of material 
being displaced into the undisturbed zone, as shown in Figure We assume 
a noncohesive Coulomb failure of the material internally, along a failure plane 
parallel to and meeting the surface of the transition zone. We additionally 
assume that the internal friction angle is also defined by Limit equilibrium 
in the x direction can then be expressed using the notation defined in Figure]^ 
as Fx = Tb + Ti COS0, where F^ = Dar, Tb = fiWo = ppVmg/{‘^■tanO) and 

Ti = pWi cos 9 = piD'^ppVrn94>o cos 0/(2 tan 0) per unit length in the y direction. 
After some rearrangement this implies that 


(7t 


pDppVrng 

2 tan0 


(1 + 00 cos^ 0). 


(5) 


This assumption of the failure surface introduces no new parameters into the 
model, and as will be shown in the following, closely predicts the measured value 
of (Tx for a large range of system parameters. In the limit where 0o —>■ 0, this 
definition reduces to that used in and m- Including this new definition of 
the threshold stress. Equation ([^, in Equation Q gives 


A best fit estimate is used to find from and is shown as the dashed 
line in Figure at a: = 0, using the measured values of I'm, 9 and 0o, which 
adequately captures the behaviour of the system past LID = 2. Before this 
limit, stress redistribution has not saturated, and is less than the predicted 
value. We find the same transition value of LjD « 2 for all cases reported here. 

The measured value of apparent friction p, = a^zlcczz inside the packing, 
shown in Figure shows that the system is far from failure inside the plug 
zone, and increasingly unstable in the transition zone. The out of plane Janssen 
coefficient, Ky = (Jyy/cjxx, shown in Figure [sp, has large variations in x, but 
is not significantly affected by the formation of the plug. The in plane Janssen 
coefficient (Figure [^), however, is strongly influenced by the 

formation of the plug, and is relatively constant with increasing L inside the 
plug zone. 
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Figure 5. Particle positions for four simulations with varying gap spacing. 
Top to bottom: D = 1, 2, 5 and 15. For all cases, (po « 0.5, I = 0.01 and the 
simulation is displayed at the time corresponding to L = 40. Particles are 
coloured by size, darker colours representing smaller particles. 



D D 


Figure 6. Descriptors of the system with varying gap spacing D. Error bars in 
each plot represent one standard deviation of the measured value. (A) Average 
solid fraction within the plug zone, Vm- (B) Slope of the pile in the transition 
zone, 9. (C) Crosses represent best fit value of Kz from measurement of the 
stress on the piston head, using Eq (|^ . Kz and Ky , represented by dots and 
triangles respectively, are calculated directly from the coarse grained granular 
packing. (D) Threshold stress ut. Dots represent the mean value of the coarse 
grained continuum field dxx a,t x = L, and crosses represent predicted values 
from Eq ([^ using measured values of r'm, 9 and po. 
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Figure 7. Descriptors of the system with varying initial packing fraction 
Error bars in each plot represent one standard deviation of the measured value. 
(A) Average solid fraction within the plug zone, v^- (B) Slope of the pile in the 
transition zone, 6. (C) Crosses represent best ht value of from measurement 
of the stress on the piston head, using Eq ([^. and Ky, represented by 
dots and triangles respectively, are calculated directly from the coarse grained 
granular packing. (D) Threshold stress ar- Dots represent the mean value of 
the coarse grained continuum field a^x x = L, and crosses represent predicted 
values from Eq using measured values of t'm, d and <j)Q. 


An underlying assumption of the Janssen stress redistribution is that when 
averaging over the width of the system (here in the y and z directions), the 
principal stress directions are parallel to the system geometry [8]. For this 
reason we measure a, the absolute value of the x component of the eigenvector 
of the major principal stress, which is shown in Figure |^. When awl the 
major principal stress points along the cc-axis, and when a « 0, the principal 
stress lies in the yz plane. For x < L we find that the major principal stress is 
collinear with the system geometry, and the Janssen stress model fits well. In a 
traditional silo problem, gravity acts parallel to the applied stress, and averaging 
across the width of a silo ensures the validity of this assumption. For this case, 
however, because the direction of gravity has broken the inherent symmetry of 
the silo problem, its validity is not ensured [T3j. We do, however, hnd that this 
assumption holds well for this and all simulations reported here. 
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3.1 Gap spacing 

As motivated in Equation Q, the gap spacing D largely controls the magnitude 
of the stresses within the system. For this reason, we here vary this spacing 
systematically from Z? = 1 to Z? = 15, while maintaining = 0-5 ± 0.05 and 
Z = 0.01 (except for the case of ZZ = 1, where ~ 0.66), as depicted for 
four values of D in Figure To make a reasonable comparison between these 
simulations, in each case the area of the piston is kept constant, such that its 
area \bWD = 100. Select measures of the behaviour of the system are shown 
in Figure Slope angles, 0, are averaged over the times corresponding to 
2D < L < lOD, whilst K and (Tt are averaged temporally over values in the 
range 2D < L < lOD, where at each time we measure spatially in the range 
L/4 < X < ZLjA:. In Figure]^ and B, we observe increasing solid fraction in 
the plug, Vrm and slope angle, 0, with increasing gap spacing, as the effect of 
the boundaries on the system decreases. For D < 8, we note that the effect of 
the 

For each simulation, the measured normal stress at the piston, a^, is fitted 
with Equation Q, and a best fit estimate of is shown as crosses, with the 
standard deviation of the error of the regression used as error bars, in Figure 
The mean and standard deviation of the measured values of and Ky 
from the continuum data between L = 2D and L = lOZZ are shown as dots and 
triangles, respectively. For D > 2 we find that both the measured values of 
and Ky are independent of D, and have mean values of Ky = 0.67 ± 0.04 and 
Kz = 0.58±0.05. Best fit estimates of K^ are also independent of D, with mean 
values of K^ = 0.56 ± 0.07. Figure]^ shows the mean of ctt also from L = 2 to 
L — 10. Triangles represent the prediction from Equation [fusing the measured 
values of Vm, 0 and (j)o. In all cases, we find the measured and fitted values of 
Kz to be in agreement, whereas the values of ar agree only with D > 3. We 
note, however, that ax depends strongly on 0, and we have as yet no means for 
predicting this quantity. The dependence of ctt on 0 is in contrast to studies 
on fold and thrust belts [3], where there is no confinement vertically above the 
material. 

3.2 Initial packing fraction 

Existing models nisi for the evolution of the system have neglected any effect 
of the initial packing fraction (pQ. To test this assumption, we here vary pQ 
from 0.1 to 0.6, while maintaining D = 5 and / = 0.01. As shown in Figure 

and B, the packing fraction inside the plug and the slope of the transition 
zone are independent of the initial packing fraction. In Figure [^, we observe 
that the measured and best fit values of Kz are in agreement for a wide range 
of 00 1 and that these values are lower than the measured values of Ky. The 
prediction of threshold stress from Equation ([^, which includes a dependence 
on 00, slightly under-predicts the threshold stress at 0o = 0.1. Nevertheless, 
both the measured and predicted values of the threshold stress in Figure [ZP 
are in agreement with observations from |3] , where a non-dimensional threshold 
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Figure 8. Descriptors of the system with varying inertial number I. Error 
bars in each plot represent one standard deviation of the measured value. (A) 
Average solid fraction within the plug zone, I'm- (B) Slope of the pile in the 
transition zone, 9. (C) Slope of the best fit value of L(s) denoted by black dots 
against prediction using incompressibilty shown as the shaded region, which 
denotes ± one standard deviation around the mean value for each case. (D) 
Crosses represent best fit value of Kz from measurement of the stress on the 
piston head, using Eq §. Kz and Ky, represented by dots and triangles 
respectively, are calculated directly from the coarse grained granular packing. 
(E) Threshold stress tr-r- Dots represent the mean value of the coarse grained 
continuum field axx at x = L, and crosses represent predicted values from Eq 
([^ using measured values oi Vm, 9 and (j)o- 
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stress of ar = 10.7 was found to reproduce the observed pattern formation 
behaviour in the quasi-static limit, at = 5, for a range of values of (jjo < 0.5. 

3.3 Inertial number 

Finally, we wish to comment on inertial effects in such a system. Towards 
this end we systematically vary the inertial number / from 10“^ to 10 while 
maintaining D = 5 and (j>o = 0.5 ± 0.05. As shown in Figure a transition 
occurs at I « 0.1, where the quasistatic behaviour begins to be dominated by 
inertial effects, and the system is fluidized. In Figure]^, we notice a jump in 
the maximum solid fraction, as the particles begin to flow and rearrange due to 
the increased piston velocity. This is accompanied by an increase in the slope 
angle 6 (Figure [^), and a decrease in the accumulation rate (Figure [^), as 
the grains begin to slip relative to the walls. With increasing piston velocity, 
the anisotropy of the system is lost, as shown in Figure |8p, and both Ky and 
Kz tend towards a mean value of 0.56 ± 0.02 for / > 1. As shown in Figure]^, 
the threshold stress, ctt, also diverges above / = 0.1 away from the theoretical 
prediction. For values of / > 0.1, we have therefore used the measured value of 
ax in the best fit estimation of shown in[8p, rather than the value predicted 
from ([^, as used in all other cases. 


4 Conclusions 

We have here described a large number of simulations of granular material which 
have been compacted in a confined geometry. For all cases, we observed that the 
stress distribution within the packing is well approximated by previous models, 
once a more rigorous definition of the threshold stress is used. This is true for 
a wide range of gap spacings, initial filling fractions and piston rates. 

In this study we have used a rough boundary condition, where macroscopic 
friction at the piston and walls is equal to the inter-particle friction. However, 
in many systems we expect the roughness at the boundaries to be lower than 
that between particles. It is unclear how this difference will affect either the 
accumulation of material near the piston head, or the stress distribution within 
the packing. 

Below a gap spacing of 3 particle diameters, the stress distribution is not 
well represented by this model. We conclude that D = 3 represents the smallest 
system size which may reasonably be described by the one dimensional Janssen 
stress model. In addition, at inertial numbers of / > 0.1, we find that there 
is significant slip at the boundary, and the threshold stress diverges from the 
model prediction. In all cases, we cannot as yet predict the slope of the free 
surface in the transition zone, but we observe that this slope approaches the 
angle of repose for large systems at low piston rates. Janssen stress coefficients 
for this system are well represented by = 0.6 ± 0.1 and Ky = 0.7 ± 0.1 for 
a wide range of system parameters. A model for the threshold stress has been 
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presented using limit equilibrium, and this holds well for systems with D > 3 
and / <0.1. 

The slope angle, 9, has been measured for different system parameters to 
lie in the range of 2°-18°. A priori^ we could only assume that this angle 
must be less than or equal to the angle of repose, which for these grains is 
approximately 20°. The wide variability of 6 is as yet unexplained, and is in 
stark contrast to the case where a top boundary does not exist, for example in 
fold and thrust belts [3]. We do note, however, that at large values of D the 
slope angle approaches the angle of repose. 

With regards to the two Janssen parameters, we can clearly distinguish the 
values of Ky and in Figure]^, Figure]^ and Figure]^, for / < 0.1. The 
reason for the difference between these two quantities may either be due to 
anisotropy in the granular packing, or due to the differing boundaries in the y 
and z directions. As the Janssen parameters are relatively insensitive to the gap 
spacing D, we conclude that this anisotropy is due to the accumulation process, 
which creates a preferential direction within the packing. This distinction is 
important when considering models which account for more complex geometries, 
such as in [3]. 
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